% This file can generate Figure 5, Figure 6, and Table 4 in the paper. 
% Also, the file is the one generating the results needed to reproduce
% Figure 8. In order to obtain the latter, this file needs to be run
% several times, one for each measure of output (option 'use_data' below);
% each time the file is being executed, it will automatically save the
% results in a .xlsx file that later will be read in the file named
% 'FIG_8_Robustness_Data_Inflation.m'. 

close all; clear; clc;

global P2 lambda PP H
%%

graph       =   1;              % = 1 to generate a plot. 

csv_file  	=   1;              % = 1 if you want to import from CSV.
                                % = 0 if you want to import from XLSX.
                                
                                
use_data    =   5;              % = 1  if inflation using Michigan Expectations
                                % Median. BASELINE    
                                % = 2  if inflation using Michigan Expectations
                                % Mean.                     
                                % = 3  if inflation using SPF CPI Median.
                                % = 4  if inflation using SPF GDP Deflator.
                                % = 5 if inflation using the Cleveland Expectations


table       =   3;              % = 1 to generate the first two columns of 
                                %     Table 4.
                                % = 2 to generate the pointwise test
                                %     (third column in Table 4).
                                % = 3 to generate the pathwise test (last 
                                %     column in Table 4)
                                 
%% Importing and Preparing the data.

if csv_file == 1
    data            =   csvread('Data_File_Baseline.csv');        
else
    data            =   xlsread('Data_File.xlsx','Baseline');
end
% Import the specific variables. 
dz          =   data(:,3);
dy          =   data(:,4);
zlb         =   data(:,5);
rec         =   data(:,6);
zlev        =   data(:,7);
ylev        =   data(:,8);
time        =   data(:,9);

if use_data == 1        % Michigan Survey of Consumer Expectations Median.
    pilev_data   =   csvread('Data_File_Inflation.csv');
    pilev        =   pilev_data(:,6);   
elseif use_data == 2    % Michigan Survey of Consumer Expectations Mean.
    pilev_data   =   csvread('Data_File_Inflation.csv');
    pilev        =   pilev_data(:,7);   
elseif use_data == 3    % SPF forecast of CPI inflation.
    pilev_data   =   csvread('Data_File_Inflation.csv');
    pilev        =   pilev_data(:,8);   
elseif use_data == 4    % SPF forecast of GDP deflator inflation.
    pilev_data   =   csvread('Data_File_Inflation.csv');
    pilev        =   pilev_data(:,10);    
elseif use_data == 5    % Cleveland's Fed expected inflation.
    pilev_data   =   csvread('Data_File_Inflation.csv');
    pilev        =   pilev_data(:,18); 
end

% Defining the sample:
st          =   148;                        % = 1 for full sample 
                                            % = 148 for 1984Q1 
                                            % = 52 for 1960Q1.                                            
                                            
%% Options for local projections.

P           =   2;      % Number of lags for endogenous variables.
H           =   9;      % Maximum number of leads.
h1          =   1;      % For splines.        
lambda      =   100;    % Value of penalty of SLP.
r           =   3;      % r-1 = order of the limit polynomial for the SLP.

%% Creating variables.

% LHS:
pilev2      =   (1/4)*pilev(st:end);
% Output:
ylev2       =   ylev(st:end);
% Productivity:
zlev2       =   zlev(st:end);
% Dummy for ZLB:
zlb2        =   zlb(st:end);
% Creating W matrix: 
ff          =   [ylev2 zlev2 pilev2];
w           =   lagmatrix(ff,1:P);  % Remove period t observations.

% Creating lagged variables interacted with period-t ZLB dummy:
[TY,YT]     =   size(w);
zlb3        =   zeros(TY,YT);
for jx = 1:YT
    zlb3(:,jx)  =   zlb2;
end
www         =   zlb3.*w; 
zlb4        =   1-zlb3;
wwww    	=   zlb4.*w;

%% Creating the IRFs.

% LHS variable:   
y                       =   pilev2(1+P:end,:); 

% Local projection.
        % Standard case.  
w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                             wwww(1+P:end,:), www(1+P:end,:)];
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); % Shock variable.   
lp                      =   locproj_var(y,x,w,h1,H,'reg'); 
if table == 3
    lp                      =   locproj_var_sur(y,x,w,h1,H,'reg'); 
end
% Extract the IRF from local projection: 
irf                     =   lp.IR(2:H+1); 
% Number of regression parameters estimated:
Xp                      =   size(lp.X,2)/H;
% Construct standard errors:
u                       =   lp.Y - lp.X*lp.theta;
[nwse_lp,vc_lp]         =   NeweyWest(u,lp.X,0,0);
    % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                             wwww(1+P:end,:), www(1+P:end,:)];
x                       =   zlb2(1+P:end).*zlev2(1+P:end,:);  
lp_zlb                  =   locproj_var(y,x,w,h1,H,'reg'); 
if table == 3
    lp_zlb                  =   locproj_var_sur(y,x,w,h1,H,'reg'); 
end
% Extract the IRF from local projection: 
irf_zlb                 =   lp_zlb.IR(2:H+1); 
% Construct standard errors:
u                       =   lp_zlb.Y - lp_zlb.X*lp_zlb.theta;
[nwse_lp_zlb,vc_lp_zlb] =   NeweyWest(u,lp_zlb.X,0,0);

% Smooth local projection.
        % Standard case.
w                       =   [zlb2(1+P:end) zlb2(1+P:end).*zlev2(1+P:end),...
                             wwww(1+P:end,:), www(1+P:end,:)];
x                       =   (1-zlb2(1+P:end)).*zlev2(1+P:end,:); 
slp                     =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
% Extract the IRF from local projection: 
irf_slp                 =   slp.IR(2:H+1);
if table == 3
    slp                     =   locproj_var_sur(y,x,w,h1,H,'smooth',r,lambda);
end
P2                      =   slp.P2;
% Construct standard errors:
u                       =   slp.Y - slp.X*slp.theta;
[nwse_slp,vc_slp]       =   NeweyWest_slp(u,slp.X,0,0);
V2                      =   slp.B*vc_slp(1:slp.K,1:slp.K)*slp.B';
se_slp                  =   sqrt(diag(V2));

        % ZLB case.
w                       =   [zlb2(1+P:end) (1-zlb2(1+P:end)).*zlev2(1+P:end),...
                             wwww(1+P:end,:), www(1+P:end,:)];
x                       =   zlb2(1+P:end).*zlev2(1+P:end,:); 
slp_zlb                 =   locproj_var(y,x,w,h1,H,'smooth',r,lambda); 
if table == 3
    slp_zlb                 =   locproj_var_sur(y,x,w,h1,H,'smooth',r,lambda);
end
% Extract the IRF from local projection: 
irf_slp_zlb             =   slp_zlb.IR(2:H+1);
% Construct standard errors:
u                       =   slp_zlb.Y - slp_zlb.X*slp_zlb.theta;
[nwse_slp_zlb,vc_slp_zlb] = NeweyWest_slp(u,slp_zlb.X,0,0);
V22                     =   slp_zlb.B*vc_slp_zlb(1:slp_zlb.K,1:slp_zlb.K)*slp_zlb.B';
se_slp_zlb              =   sqrt(diag(V22));

%% Saving the IRFs and the corresponding standard errors. 

t                       =   0:H-1;
p                       =   0.9;
bb                      =   tinv(p,TY-P);

if use_data == 1
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Inflation.xlsx',save_LP,'UM MEDIAN','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_LP_ZLB,'UM MEDIAN ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Inflation.xlsx',save_SLP,'UM MEDIAN','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_SLP_ZLB,'UM MEDIAN ZLB','E2');
elseif use_data == 2
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Inflation.xlsx',save_LP,'UM MEAN','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_LP_ZLB,'UM MEAN ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Inflation.xlsx',save_SLP,'UM MEAN','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_SLP_ZLB,'UM MEAN ZLB','E2');
elseif use_data == 3
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Inflation.xlsx',save_LP,'SPF CPI MEDIAN','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_LP_ZLB,'SPF CPI MEDIAN ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Inflation.xlsx',save_SLP,'SPF CPI MEDIAN','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_SLP_ZLB,'SPF CPI MEDIAN ZLB','E2');  
elseif use_data == 4 
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Inflation.xlsx',save_LP,'SPF DEFLATOR MEDIAN','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_LP_ZLB,'SPF DEFLATOR MEDIAN ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Inflation.xlsx',save_SLP,'SPF DEFLATOR MEDIAN','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_SLP_ZLB,'SPF DEFLATOR MEDIAN ZLB','E2');
elseif use_data == 5 
    save_LP     = [irf (irf-bb*nwse_lp(1:H)) (irf+bb*nwse_lp(1:H))];
    xlswrite('Inflation.xlsx',save_LP,'CLEVELAND','B2');
    save_LP_ZLB = [irf_zlb (irf_zlb-bb*nwse_lp_zlb(1:H)) (irf_zlb+bb*nwse_lp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_LP_ZLB,'CLEVELAND ZLB','B2');
    save_SLP     = [irf_slp (irf_slp-bb*se_slp(1:H)) (irf_slp+bb*se_slp(1:H))];
    xlswrite('Inflation.xlsx',save_SLP,'CLEVELAND','E2');
    save_SLP_ZLB = [irf_slp_zlb (irf_slp_zlb-bb*se_slp_zlb(1:H)) (irf_slp_zlb+bb*se_slp_zlb(1:H))];
    xlswrite('Inflation.xlsx',save_SLP_ZLB,'CLEVELAND ZLB','E2');
end
%% Creating graph. 

if graph == 1
    figure
    plot(t,irf,'x-b',t,irf_zlb,'*-r','Linewidth',1.5);
    h = legend('No ZLB','ZLB','Location','NorthWest');
    set(h,'Interpreter','latex');
    hold on
    shadedplot(t,(irf+bb*nwse_lp(1:H))',(irf-bb*nwse_lp(1:H))',[0.75,0.75,0.75],'w');
    hold on
    plot(t,irf,'x-b',t,irf_zlb,'*-r',t,irf_zlb+bb*nwse_lp_zlb(1:H),'*--r',t,irf_zlb-bb*nwse_lp_zlb(1:H),'*--r','Linewidth',1.5);
    grid on
    title('Expected Inflation to Productivity','Interpreter','latex','fontSize',14);
    xlabel('Horizon','Interpreter','latex','fontSize',12);
    ylabel('Percent','Interpreter','latex','fontSize',12);
        
    figure
    plot(t,irf_slp,'x-b',t,irf_slp_zlb,'*-r','Linewidth',1.5);
    h = legend('No ZLB','ZLB','Location','NorthWest');
    set(h,'Interpreter','latex');
    hold on
    shadedplot(t,(irf_slp+bb*se_slp(1:H))',(irf_slp-bb*se_slp(1:H))',[0.75,0.75,0.75],'w');
    hold on
    plot(t,irf_slp,'x-b',t,irf_slp_zlb,'*-r',t,irf_slp_zlb+bb*se_slp_zlb(1:H),'*--r',t,irf_slp_zlb-bb*se_slp_zlb(1:H),'*--r','Linewidth',1.5);
    %axis([0 8 min(lb_SE)-0.2 max(up_SE)+0.2 ]);
    grid on
    title('Expected Inflation to Productivity','Interpreter','latex','fontSize',14);
    xlabel('Horizon','Interpreter','latex','fontSize',12);
    ylabel('Percent','Interpreter','latex','fontSize',12);
end
%% Creating Table 4. 
 
if table == 1
    % Point estimates and corresponding standard errors. 
    
    Beta_h              =   [irf_slp se_slp];
    Beta_h_z            =   [irf_slp_zlb se_slp_zlb];

    Headers             =   {'Beta_h', 'Std_Errors'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'};   
                         
    Table_3_Beta_h      =   array2table(Beta_h,'RowNames',Rowish,'VariableNames',Headers)
    
    Headers_z           =   {'Beta_h_z', 'Std_Errors'};
    Table_3_Beta_h      =   array2table(Beta_h_z,'RowNames',Rowish,'VariableNames',Headers_z)   

elseif table == 2
    % Pointwise test.
    
    Fslp                =   zeros(H,2);
    Flp                 =   zeros(H,2);
    chislp              =   zeros(H,2);
    
    for jx = 1:H
        Fslp(jx,1)          =   ((irf_slp(jx) - irf_slp_zlb(jx))^2/(V2(jx,jx) + V22(jx,jx)));
        Flp(jx,1)           =   ((irf(jx) - irf_zlb(jx))^2/(vc_lp(jx,jx) + vc_lp_zlb(jx,jx)));
    end
    
    chislp(:,1)         =   Fslp(:,1);
    
    for jx = 1:H
        Fslp(jx,2)          =   fcdf(Fslp(jx),1,TY-P-jx-Xp,'upper');
        Flp(jx,2)           =   fcdf(Flp(jx),1,TY-P-jx-Xp,'upper');
        chislp(jx,2)        =   chi2cdf(chislp(jx,1),1,'upper');
    end
    
    COEFF_lp            =   [irf irf_zlb];
    SE_lp               =   [nwse_lp(1:H) nwse_lp_zlb(1:H)];
    
    Headers             =   {'Test', 'P_value'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'};   
    Table_3_Pointwise   =   array2table(chislp,'RowNames',Rowish,'VariableNames',Headers)
    
elseif table == 3
    % Pathwise test.
        % Construct combined variance-covariance matrix
    Vtog                =   zeros(2*H,2*H);
    Vtog(1:H,1:H)       =   V2;
    Vtog(H+1:end,H+1:end) = V22;
    R                   =   zeros(H,2*H); 
    R(:,1:H)            =   eye(H,H);
    R(:,H+1:end)        =   -eye(H,H);
    rr                  =   zeros(H,1);
    betahat             =   [irf_slp;irf_slp_zlb];
    
    for jx = 1:H
        Fslp_sur(jx,1)      =   ((irf_slp(jx) - irf_slp_zlb(jx))^2/(V2(jx,jx) + V22(jx,jx)));
        Flp_sur(jx,1)       =   ((irf(jx) - irf_zlb(jx))^2/(vc_lp(jx,jx) + vc_lp_zlb(jx,jx)));
    end
    
    for jx = 1:H
        Fslp_sur(jx,2)      =   fcdf(Fslp_sur(jx),1,TY-P-jx-Xp,'upper');
        Flp_sur(jx,2)       =   fcdf(Flp_sur(jx),1,TY-P-jx-Xp,'upper');
    end
    
    Fpath_sur            =   zeros(H,2);
    chipath_sur          =   zeros(H,2);
    for jx = 1:H
        RR                  =   R(1:jx,:);
        Fpath_sur(jx,1)     =   (RR*betahat)'*(RR*Vtog*RR')^(-1)*(RR*betahat);
    end
    chipath_sur(:,1)    = 	Fpath_sur(:,1);
    
    for jx = 1:H
        Fpath_sur(jx,2)     =   fcdf(Fpath_sur(jx,1),jx,TY-P-H-Xp,'upper');
        chipath_sur(jx,2)   =   chi2cdf(chipath_sur(jx,1),jx,'upper');
    end
            
    Headers             =   {'Test', 'P_value'};
    Rowish              =   {'h_equal_0','h_equal_1','h_equal_2','h_equal_3',...
                             'h_equal_4', 'h_equal_5', 'h_equal_6','h_equal_7',...
                             'h_equal_8'};   
    ColumnResults 		= 	[chipath_sur(:,1) Fpath_sur(:,2)];                        
    Table_3_Pathwise    =   array2table(ColumnResults,'RowNames',Rowish,'VariableNames',Headers)

end
    

    